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Abstract Tethering methods allow us to perform Monte Carlo simulations in ensem- 
bles with conserved quantities. Specifically, one couples a reservoir to the physical 
magnitude of interest, and studies the statistical ensemble where the total magnitude 
ri ■ (system+reservoir) is conserved. The reservoir is actually integrated out, which leaves 

O ' us with a fluctuation-dissipation formalism that allows us to recover the appropriate 

^5 . Helmholtz effective potential with great accuracy. These methods are demonstrating a 

P^ ' remarkable flexibility. In fact, we illustrate two very different applications: hard spheres 

crystallization and the phase transition of the diluted antiferromagnet in a field (the 
Cu ' physical realization of the random field Ising model). The tethered approach holds 

C/3 , the promise to transform cartoon drawings of corrugated free-energy landscapes into 

real computations. Besides, it reduces the algorithmic dynamic slowing-down, probably 
because the conservation law holds non-locally. 



Keywords Monte Carlo methods, barriers, effective potential 
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O ■ 1 Introduction 

Monte Carlo simulation is one among the handful of general methods that physicists 

^ ' can use to explore strongly coupled problems, far from the perturbative regime [I] . Fur- 

V^ I thermore, the Monte Carlo method has become itself an object of active investigations. 

0^ ■ The reason is twofold: since it is one of our most cherished tools we wish to sharpen it. 

>— -^ ' But, it is also true that the dynamic bottlenecks found during the simulation resemble 

closely the dynamic arrests that one may find in Nature. 
P^ ' Here, we discuss a strategy to address a fairly general problem in Monte Carlo 

>— ' , simulations: free-energy barriers. Indeed, the time that a standard, canonical simulation 

gets trapped in a local minimum grows exponentially with the free-energy barrier to 

be surmounted. Let us recall some important instances of this generic problem (the list 

is far from exhaustive): 

" T— I . 
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In every first-order transition tlie system spontaneously segregates in a spatially 
heterogeneous mixture of the two phases (think, e.g., of ice and water at the melting 
point). The two coexisting phases are separated by an interface. As a consequence, 
in order to produce a significant change, the system must build an interface of 
linear dimensions comparable to those of the simulation box, L. The corresponding 
free-energy cost scales with the system's cross-section L , where D is the space 
dimension [21. As a consequence, r, the simulation characteristic time grows as 
r oc expfX"!/ ~ ], where X" is a surface tension. This disaster is named exponential 
dynamic slowing down. 

Studies of crystallization, the first-order phase transition encountered upon cool- 
ing or compressing a liquid, are hampered by problems worse than exponential 
dynamic slowdown. Even for such a simple model liquid as hard spheres, there is 
an impressively large number of free-energy minima where the simulation may get 
stuck. For instance, although for simple model liquids the equilibrium crystal is 
face-centered cubic, the system might be prone to nucleate a body-centered cubic 
phase [3]. Even if a crystal of correct symmetry is formed, it may have a large 
number of defects (point defects, such as vacancies, or non-local defects such as 
dislocations). Furthermore, an amorphous solid (a glass) may appear [4]. 
Even if the phase transition is continuous, large barriers might be present at the 
critical temperature, as it is the case for the random field Ising model [5], where 
AF Qc L <^ L ~ . Note that, for this problem, both the appropriate order param- 
eter and the microscopic ground state configuration ^ are known. However, the 
simulation gets trapped in local minima with escape times logr ~ L . An added, 
major difficulty is the need to simulate a huge number of samples to obtain disorder 
averages truly representative of the system's behavior. 

In spin glasses [7] not even the appropriate order parameter is known. To detect 
the phase transition one must use real replicas (clones of the system, evolving un- 
der the same Hamiltonian but with uncorrelated thermal noise). There is a fairly 
large number of local minima where the simulation may trap (at least for finite sys- 
tems 8j). Some specific Monte Carlo methods have been devised for this problem, 
such as the exchange Monte Carlo method [51llU| (also known as parallel temper- 
ing). Furthermore, specific hardware has been used to study these systems [1111121 
fT3l[T4l[l 5, 16, 17, 18.19.20.21 . Unfortunately, exchange Monte Carlo suffers a strong 
(probably exponential) dynamic slowing down below the critical temperature [2U| . 



We have some effective and general methods to cope with the first type of problem, 
namely first-order phase transitions where the high and the low-temperature configura- 
tions are clearly identified, and easily findable. Multicanonical [22] or Wang-Landau [23] 
simulations feature a generalized statistical ensemble. The system performs a random 
walk in the energy space, back and forth from the energy of the ordered phase to that 
of the disordered phase. At the price of optimizing a number of parameters propor- 
tional to the system size, the free-energy barrier separating these two phases is easily 
overcome in small systems. However, the random-walk strategy can only delay to larger 
system sizes the advent of exponential dynamic slow down. In fact, for large enough 
systems, surface-energy effects induce geometric transitions in the energy gap between 
the ordered and the disordered phase [24.25,26,27,28 . The energy random walk needs 
a huge amount of simulation time in order to tunnel through the geometric transitions, 
which results in exponential dynamic slowing down |29j . 



A difFerent approach stems from Lustig's microcanonical Monte Carlo |30] . Given its 
microcanonical nature, one may perform independent simulations at different energies 
located in the gap between the ordered and the disordered phase. In this way, one 
avoids the tunneling through the geometric transitions, which allows the equilibration 
of larger systems. Using a fluctuation-dissipation formalism one can recover the entropy 
density, from which a precision study of the phase transition follows [31] |J 

However, in general applications, the gap between coexisting phases is not an en- 
ergy gap. Rather, some reaction coordinate (such as an order parameter) labels the 
different regions of the phase space that we wish to explore. For instance, liquid or 
crystalline configurations in hard spheres all have the same energy. Just as entropy 
is the natural thermodynamic potential when the reaction coordinate is the internal 
energy, in the general case one wishes to study the Helmholtz effective potential associ- 
ated to the appropriate reaction coordinate. This is the case, for instance, in the study 
of crystallization kinetics in supercooled liquids ._3i,33j, where the reaction coordinate 
is a bond-orientational crystalline order parameter [35 ■ A variation of the random- 
walk strategy to reconstruct the Helmholtz potential, named umbrella sampling [35j . 
is popular in this context. 

A recent alternative is Tethered Monte Carlo 36 , which allows one to reconstruct 
the Helmholtz effective potential without random walks in the space of the effective co- 
ordinate. The method is a generalization of Lustig's microcanonical Monte Carlo [30| . 
One simulates an ensemble where the reaction coordinate is (almost) constrained to 
take any desired value. The associated fluctuation-dissipation formalism permits the ac- 
curate reconstruction of the Helmholtz potential and (if it is so wished) of the canonical 
expectation values. Nevertheless, one is not restricted to a mere speed-up of a canoni- 
cal simulation. Very interesting new information can be unearthed from this different 
viewpoint. 

Furthermore, for a complex enough system, a single reaction coordinate would not 
suffice. Wang-Landau, or umbrella sampling simulations are rather cumbersome if one 
needs to fine-tune parameters for a two-dimensional random walk. On the other hand, 
as we shall show below, a tethered approach to the problem is not necessarily tougher 
than the one-dimensional case. 

Our aim here is to describe Tethered Monte Carlo, with an emphasis on its sim- 
plicity and generality. We illustrate it with its application to two challenging problems: 
hard-spheres crystallization and the diluted antiferromagnet in an external field (the 
much easier problem of the ferromagnetic Ising model was considered in Refs. [36II37| ). 
We make emphasis on the algorithmic aspects of the computation. For a discussion of 
the physical results, the reader is referred to Refs. [38,39,40,41 . It turns out that in 
both problems one needs to compute a Helmholtz potential depending on two order 
parameters. However, the crystallization transition is of the first order, while that of 
the diluted antiferromagnet is a continuous one. Hence the strategy followed in each 
case differs. 

Let us mention as well a tethered computation of the Helmholtz effective potential 
for a simple model of glass- forming liquid [42]. A most remarkable feature is that, for 
this system, an appropriate reaction coordinate is unknown. However, the problem can 
be bypassed considering real replicas. 



^ It was emphasized by W. Janke that the entropy is a more natural thermodynamic poten- 
tial than the free energy for the analysis of first-order transitions |32| . 



The layout of the rest of this paper is as follows. Section [2] describes the general 
set-up for a tethered simulation using a simple Metropolis update, in the familiar con- 
text of the Ising model. We note en passant that, if a Kasteleyn-Fortuin decomposition 
is known for the canonical version of the problem at hand, a cluster update can be 
implemented for the Tethered Monte Carlo [37] • A detailed derivation of the tethered 
formalism is given in Sect.[3l In Sect. [4] we consider the first problem where two con- 
served quantities are needed, the diluted antiferromagnet in a field. Sect. [S] describes 
the application of the tethered formalism to hard-spheres crystallization. We give our 
conclusions in Sect. (6] 



2 Tethered Monte Carlo, in a nutshell 

In this section we give a brief overview of the Tethered Monte Carlo (TMC) method, 
including a complete recipe for its implementation in a typical problem. This is as sim- 
ple as performing several independent ordinary MC simulations for different values of 
some relevant parameter and then averaging them with an integral over this parameter. 
We shall give the complete derivations and the detailed construction of the tethered 
ensemble in Section [31 

We are interested in the scenario of a system whose phase space includes several 
coexisting states, separated by free-energy barriers. The first step in a TMC study is 
identifying the reaction coordinate x that labels the different relevant phases. This can 
be (but is not limited to) an order parameter. In the remainder of this section we shall 
consider a ferromagnetic setting, so the reaction coordinate will be the magnetization 
density m. 

The goal of a TMC computation is, then, constructing the Helmholtz potential 
associated to m, 0]^[l3,m), which will give us all the information about the system. 
This involves working in a new statistical ensemble tailored to the problem at hand, 
generated from the usual canonical ensemble by Legendre transformation: 

ZN[P,h) = e^^"('^^'^) = /dm ^mhm^O^{P,ra)\ _ ^^^ 



where Zj^{P,h) is the canonical partition function, F]\[{j3,h) the Gibbs free-energy 
density and A^ is the number of degrees of freedom. 

Since in a lattice system the magnetization is discrete, we actually couple it to a 
Gaussian bath to generate a smooth parameter, called m. The effects of this bath are 
integrated out in the formalism. 

In order to implement this construction as a workable Monte Carlo method we need 
to address two different problems: 

— We need to know how to simulate at fixed m. 

— We need to reconstruct Qj^{P,'m) from simulations at fixed m, and afterwards, to 
recover canonical expectation values from Eq. ([T]) to any desired accuracy. 

We shall explain separately how to solve each of the two problems, in the following two 
paragraphs. 



2.1 Metropolis simulations in the tethered ensemble 

Let us denote the reaction coordinate by m (for the sake of concreteness let us think 
of the magnetization density for an Ising model) . The dynamic degrees of freedom are 
{si\ (they could be spins, or maybe atomic positions). Therefore m is a dynamical 
function (i.e. a function of the {sj}). We wish to simulate at fixed rh (rh is a parameter 
closely related to the average value of m) . 

The canonical weight at inverse temperature j3 would be exp[— /3t/] where U is 
the interaction energy. Instead, the tethered weight is (see Ref. ^36, and Sect. |3] for a 
derivation) 

a;jv(/?, m; {s,}) = ^-PU+N(ra-n^) ^^^ _ ^)(Ar-2)/2^(^ _ ^) (2) 

The Heaviside step function 9(rh — m) imposes the constraint that rh > m({sj}). 

The tethered simulations with weight ^ are exactly like a standard canonical 
Monte Carlo in every way (and satisfy detailed balance, etc.). For instance, in an Ising 
model setting, the common Metropolis algorithm l43l is 

1. Select a spin Sj. 

2. The proposed change is flipping the spin, Si — >■ — s^.U 

3. The change is accepted with probabilitj|j 

P{si — 5> —Si) = min{l,ajncw/<.<Jold}- (3) 

4. Select a new spin Sj and repeat the process. 

We remark that the above outlined algorithm produces a Markov chain entirely 
analogous to that of a standard, canonical Metropolis simulation. As such it has all the 
requisite properties of a Monte Carlo simulation (mainly reversibility and ergodicity). 
Tethered mean values can be computed as the time average along the simulation of 
the corresponding dynamical functions (such as internal energy, magnetization den- 
sity, etc.). Statistical errors and autocorrelation times can be computed with standard 
techniques [44ll45] . 

The actual magnetization density is constrained (tethered) in this simulation, but 
it has some leeway (the Gaussian bath can absorb small variations in m, see Sect.[31). 
In fact, its fluctuations are crucial to compute an important dynamic function, whose 
introduction would seem completely unmotivated from a canonical-ensemble point of 
view: 

-._ 1 aioga;jv(/3,m;{s,}) _ JV - 2 

iV dm 2N[m-m{{si})\' ^' 

One of the main goals of a tethered simulation is the accurate computation of the 
expectation value {h)rh- 



^ In an atomistic simulation, one would try to displace a particle, or maybe to change the 
volume of the simulation box. 

^ In general, the Metropolis probability has to take into account both the weight and the 
probability of proposing this particular change U. However, for this simple problem the latter 
is symmetric, so it quotients out. 



The case where one wishes to consider two reaction coordinates mi and m,2 is 
completely analogous: 

UMW,mi,m2;{s,}) = ^-0U+N{m^-rh^)+N{m2-rh2) ^ (5) 

X (mi — mi)^ ^ '' (m2 — m,2) 9{mi — mi)9{iri2 — 1712) , 

f ^ 1 d log 0JNil3,mm2;{si}) ^ A^ - 2 

^~ iV 9mi 2Af[mi -mi({s,})] ' ^ '' 

? ^ 1 91ogajjv(/3,m,m2;{sJ) ^ A^ - 2 

2- iV 9m2 2iV[m2-m2({si})l • ^' 

We note as well that, in the context of supercooled liquids crystallization [3], a 
tethered weight slightly different from Eq. © has been in use (however, the importance 
of the dynamic function b was apparently not recognized) : 

^N[P,rn-{si}) = e ^ "• ' ' , (8) 

b = a{rh — m) . (9) 

Here, a is a tunable parameter that can be handy to control the deviations of m from 
m [31]. Note as well that the weight ([SJ does not impose rh > m. 

For the Ising model, a Metropolis Tethered Monte Carlo simulation reconstructs 
the crucial tethered magnetic field 6 without critical slowing down (see 36 for a bench- 
marking study). This may be considered surprising for what is a local update algorithm, 
but notice that the constraint on m is imposed globally. Non-magnetic observables, 
such as the energy, do not enjoy this non-local information and hence show a typical 
2 « 2 critical slowing down (although the correlation times are low enough to permit 
equilibration for very large systems, ,36j). 

Let us stress that the above outlined update algorithm is by no means the only one 
possible. For instance, the Fortuin-Kasteleyn construction [461147) can be performed 
just as easily in the tethered ensemble, so we can consider tethered simulations with 
cluster update methods [4811491150] . This is demonstrated in [37], where the tethered 
version of the Swendsen-Wang algorithm is shown to have the same critical slowing 
down as the canonical one for the D = 3 Ising model {z ~ 0.47). This is an example 
that the use of the tethered formalism implies no constraints on the choice of Monte 
Carlo algorithm, nor hinders it in the case of an optimized method. 



2.2 Reconstructing the Helmholtz effective potential from simulations at fixed m 

The steps in a TMC simulation are, then, (see also Figure [T]) 

1. Identify the range of rh that covers the relevant region of phase space. Select A'^^ 
points rfii, evenly spaced along this region. 

2. For each rhi perform a Monte Carlo simulation where the smooth reaction param- 
eter m will be fixed at m = rhi. 

3. We now have all the relevant physical observables as discretized functions of rh. 
We denote these tethered averages at fixed rh by {O)^. 

4. The average values in the canonical ensemble, denoted by (O), can be recovered 
with a simple integration 



L 



{0)= / dmp(m)(0)^. (10) 




Fig. 1 (Color online) Computation of the Helmholtz potential i7jv and the canonical expec- 
tation values from tethered averages in a -D = 3, L = 64 ferromagnetic Ising model at the 
critical temperature. The upper panel shows the tethered magnetic field (b)mi with an inset 
zooming in on the region between its two zeros. The statistical errors cannot be seen at this 
scale (except for the leftmost points in the inset). The integral of this quantity is the Helmholtz 
potential i?jv("i)- The middle panel shows the p{m) = exp[— N f2 1^ {rh)] in a linear (right axis) 
and in a logarithmic scale (left axis) . Finally, the bottom panel shows the tethered expectation 
values of the energy density u. Their integral over the whole rh range, weighted with p{m), 
gives the canonical expectation value (u) = —0.3322894(36) (horizontal line). See Ref. 1371 for 
further details on these simulations. 



In this equation the probability density p(m) is 

p(m) = e~^^^"*^™\ r2Af(m) = r2Af(m,nin)+ / dm' (6)^.. (11) 

The tethered field {&)m was defined in Eq. Q. The integration constant fiAr(rnmin) 
is chosen so that the probability is normalized. 

5. If we are interested in canonical averages in the presence of an external magnetic 
field h, we do not have to recompute the (O)^, only O^^. This is as simple as 
shifting the tethered magnetic field: (6)m ~^ {b)m ~ Ph. 

6. In order to improve the precision and avoid systematic errors, we can run additional 
simulations in the region where p(m) is largest. 

The whole process is illustrated in Figure [U where we compute the energy density 
at the critical temperature in an L = 64 lattice of the D = 3 Ising model. Notice 
that the tethered expectation values {u}^ vary in about 10% in our m range, but 
the computation of the effective potential is so precise that the averaged value for the 
energy, {u) = —0.3322894(36), has a relative error of only ~ 10~ . 



Table 1 Energy density of the D = 3 ferromagnetic Ising model computed with the Tethered 
Monte Carlo method, showing that the reconstruction of canonical averages can be performed 
with great accuracy. The second column shows the number of points in the rh grid and the 
third the number of Monte Carlo sweeps taken on each (we use a cluster update scheme). 
See |37| for further details on these simulations. 

L Nrh MCS -{u) 

IE 91 W 0.344905(35) 

32 91 10^ 0.335730(26) 

64 109 10'^ 0.3322894(36) 

128 50 Kf 0.3309831(15) 



This is the general TMC algorithm for the computation of canonical averages from 
the Helmholtz potential. As we shall see in some of the applications, sometimes the 
integration over all phase space in step 4 is not needed and one can use the ensemble 
equivalence property to recover the (O) from the (O)^ through saddle-point equations, 
remember Eq. IT}. In other words, the tethered averages can be physically meaningful 
by themselves. We note as well that our crystallization study in Sect. [5] is built entirely 
over the effective-potential, one never uses the p(rh). 

As will be shown in the next section, the reconstruction of canonical averages from 
the combination of tethered averages does not involve any approximation. We can 
achieve any desired accuracy, provided we use a sufficiently dense grid in rh (to control 
systematic errors) and simulate each point for a sufficiently long time (to reduce statis- 
tical ones). Table[T]shows the kind of precisions that we can achieve. One could initially 
think that the computation of the exponential in plm) — exp[—Nn]\[{rri)] would pro- 
duce unstable or imprecise results for large system sizes. Instead, the combination of 
self-averaging and no critical slowing down makes the numerical precision grow with 
N. 



3 The tethered ensemble 

In this section we construct the statistical ensemble that supports the TMC method. We 
shall consider the case of the D-dimensional ferromagnetic Ising model. Reproducing 
the steps of this derivation for any other system is straightforward. We include details 
on how to introduce several tethered variables (subsection 13. l|l . Subsection l3.2l discusses 
the relationship between tethered and canonical expectation values from the point of 
view of the ensemble equivalence property. 

The Ising model is characterized by the following partition function 



Z = ^ exp /? ^ 



{s.} 



{x,y) 



(12) 



where the angle brackets indicate that the sum is restricted to first neighbors and the 
spins are Sx = ±1. In this section we will be working at fixed /3. Hence, to lighten the 
expressions we shall drop the explicit /3 dependencies. 
The energy and magnetization of this system are 

U^Nu = -^SxSy, M = Nm^'^Sx, (13) 



(A'^ = L is the number of spins in the system). 

The canonical average of a generic observable O is 



(O) = i ^ Oas^}) e- 



pu 



(14) 



U^} 



Since this is a ferromagnetic system, we may be interested in considering the average 
value of O conditioned to different magnetization regions. The naive way of doing this 
would be 



(0|m) = 



(0^(^-Ex^x/JV)) 



The canonical average could then be recovered by a weighted average of the {0\m) 



(15) 



(O) =^(0|m)pi(m), pi(m) = (slm-J^^o^/N 



(16) 



In the thermodynamical limit, the reaction coordinate m varies continuously from —1 
to 1. For a finite system, however, there are only A'^ + 1 possible values of m, so pi{m) 
is a comb-like function. 

We want to construct a statistical ensemble where a smooth Helmholtz potential 
can be defined in finite lattices. In order to do this, the first step is extending the 
configuration space with a bath R oi N Gaussian demons and defining the smooth 
magnetization M, 

M ^ Nm = M + R. (17) 

These demons can be introduced either linearly or quadraticallylJ 






i 



(18) 
(19) 



Now our partition function is 



Z = 



P2{r) = 



/CX) / ^ 
n 



dr/i 



^exp -pU-J2vf/2 



N 

n 



dVi 



exp 



■E''?/^ 



S{r- R/N) 



(20) 



(21) 



Notice that the demons are statistically independent from the spins. The convolution 
of pi(m) and P2ij') then gives the probability density function for m. 



p{rh) = I dmdr pi{m)p2{r)5{m — m — r). 



(22) 



So p(m) is essentially a smooth version of pi(m). We can do an analogous construction 
for off-lattice systems. In this case, the demons would increase the dispersion in the 



* In fact, they need not even be Gaussian variables, although we shall not consider other 
distributions here. 
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already continuous pi{x). We can control these fluctuations by considering a tunable 
number of demons, see (I37|l and Section [5] 

We can now rewrite the canonical average (O) as 



{O) = J dm (O)^ p{m), 



(23) 



where the tethered average {O)^^ is 

(0),j, = -— -/ lld7^,J20i{s.})e~^"'-^^^^/^5{m-m-R/N). (24) 
Zp{m) J_^ A± ^ 

This way, each (O)^ will have contributions from all the (Ojm), with a weight that 
will depend on the distance between m and rh. This generates smooth (O)™ and p{rh), 
so we can now define the Helmholtz effective potential in our finite lattice as 

p(m) = e-^^^"(") (25) 

It may seem that with our extended configuration space, we have sacrificed too much 
in order to get smooth functions, but we can actually integrate the demons out with 
the Dirac deltas: 

^-Nf2.(nr)^^J2oJN{m,{s^}), (26) 

where the weight iU]\j{m,{sa;}) is, for our linear and quadratic demons (neglecting 
irrelevant constant factors), 

^^i?H^, {s.}) « e-^^+^'^-*^(m - m)(~-2)/2^(^ _ ^^ ^27) 

c.(,^)(m,{..})oce-^^-(*^-^^)^/(2^). (28) 

In general, u)js[(rh, {sx}) will be of the following form, 

a;jv('", {sx}) oc e^^'^7(m,m), (29) 

that is, the thermal component typical of a canonical simulation and a magnetic factor 
that depends only on m and rh. The weight uij^ defines a new statistical ensemble, so 
we can now rewrite the tethered averages as 

A particularly important tethered average is the m-derivative of the effective potential, 
the tethered magnetic field {&)„, 

For our two kinds of demons this is 

S(Q)^l-^/^~^/^, (32) 



h(^) 



rh — m. (33) 
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Definition pi[) . together with the condition that p(m) — e^ iv(ni) j^g normalized, 
allows us to reconstruct the Helmholtz potential from tethered averages alone. The 
canonical averages can then be recovered through Eq. (|23p . 

The tethered magnetic field is essentially a measure of the fiuctuations in m, which 
illustrates the dual roles the magnetization and magnetic field play in the canonical 
and tethered formalism. This is further demonstrated by the tethered fluctuation- 
dissipation formula 

^^ = (12) , + ^i(Ob),n - {0),n{b)n^]■ (34) 

Notice that the values of m where (b)^ = will define maxima and minima of the 
effective potential and, hence, of the probability pijh). 

The tethered ensemble is, thus, constructed via a Legendre transformation of the 
canonical one, that exchanges the roles of the magnetic field and the magnetization. 
In fact, we can write the canonical partition for a given value of the applied magnetic 
field in the following way 

ZW =e^^"('') = /dm e-^I^^W-^"™]. (35) 



One interesting consequence is that computing canonical averages for non-zero mag- 
netic fields does not imply changing (recomputing) the tethered averages, only their 
weights, 

{0)(h) = ^ |dm (0)a e-^[^^"('')-^''"l, (36) 

in other words, we just have to shift the tethered magnetic field, (6)m — >■ (&)r?i ^ Ph. 

A final general point concerns the introduction of the Gaussian demons. Our def- 
initions psp and p9p and subsequent construction make the assumption that there 
are as many demons as spins. While this choice seems natural, it is by no means a 
necessity. In fact, in order to control the fluctuations in the reaction coordinate, it may 
in occasion be interesting to consider a variable number of demons (see Section [S|: 

aN 

Ra = -y^Vi, M = Nm = M + Ra, (37) 

i=l 

where a is the tunable parameter. All the computations in this section can be easily 
reconstructed for this more general case. As an example, the tethered magnetic field 
with qA'' linearly added demons is 

ll-^^ =a{m-m). (38) 



3.1 Several tethered variables 

Throughout this section we have considered an ensemble with only one tethered quan- 
tity. However, as we shall see in the following sections, it is often appropriate to consider 
several reaction coordinates at the same time. The construction of the tethered ensem- 
ble in such a study presents no difficulties. We start by coupling reaction coordinates 
Xi, i = l..n, with N demons each, 

Xi = Nxi ^ Xi + Ri, ..., Xn = NXn ^ Xn + Rn- (39) 
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We then follow the same steps of the previous section, with the consequence that the 
tethered magnetic field is now a conservative field computed from an n-dimensional 
potential ^^^{x) 

vnN = {di,nN,...,di„nN)^B (40) 

B = ((6i)£,...,(6„)i), (41) 

and each of the 6j is of the same form as in the case with only one tethered variable. 
Similarly, the tethered weight of Eq. (I29|l is now 

^ (&, {sx}) = e"'^^ -yixi, xi)'y{x2, X2) ■ ■■'y{xn,Xn). (42) 



3.2 Ensemble equivalence and spontaneous symmetry breaking 

Throughout this section we have implicitly assumed that the final goal of a Tethered 
Monte Carlo computation is eventually to reconstruct the canonical averages. As we 
shall see in some of the applications, however, this is not always the case. One example 
is the computation of the hyperscaling violations exponent of the RFIM in [3S], which 
is made directly from i^^v- 

Still, most of the time the averages in the canonical ensemble are the ones with 
a more direct physical interpretation (fixed temperature, fixed applied field, etc.). In 
principle, their computation implies reconstructing the whole effective potential fi]\[ 
and using it to integrate over the whole coordinate space, as in (|23p or, more gen- 
erally, (|36|l . Sometimes, however, the connection between the tethered and canonical 
ensembles is easier to make. Let us return to our ferromagnetic example, with a single 
tethered quantity m. Recalling the expression of the canonical partition function in 
terms of i?7v for finite h, Eq. (|35|l . a steepest-descent tells us that, for large A'', the 
integral is dominated by one saddle point such that 

d^,[f2N{m)~l3hm]^0 ^ {b),n = Ph. (43) 

Clearly, this saddle point (actually, the global minimum of the effective potential) 
rapidly grows in importance with the system size N, to the point that we can write 

hm (0)(ft)= hm (0),j,(;,). (44) 

That is, in the thermodynamical limit we can identify the canonical average for a 
given applied magnetic field h with the tethered average computed at the saddle point 
defined by h (which is nothing more than the point where the tethered magnetic field 
coincides with the applied magnetic field). 

This ensemble equivalence property would be little more than a curiosity if it were 
not for the fact that the convergence is actually extremely fast (see [3B] for a study). 
Therefore, in many practical applications the equivalence (I44p can be made even for 
finite lattices (see Section 2] below for an example). 

The saddle-point approach can be applied to systems with spontaneous symmetry 
breaking. In the typical analysis, one has to perform first a large-A limit and then 
a small-field one (this is troublesome for numerical work, where one usually has to 
consider non-analytical observables such as |m|). In the tethered ensemble we can 
implement this double limit in an elegant way by considering a restricted range in the 
reaction coordinate from the outset (see [5S] for a straightforward example in the Ising 
model and Section [4.41 for a more subtle one). 
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4 The diluted antiferromagnet in a magnetic field 

In this section we consider tlie application of tiie tethered formalism to the diluted an- 
tiferromagnet in an external field (DAFF), one of the simplest physical systems show- 
ing cooperative behavior in the presence of quenched disorder. A good experimental 
introduction to this system can be found in [5T]. It was shown in |52ll53j that the 
DAFF belongs to the same universality class as the random field Ising model (RFIM), 
a system more amenable to analytical treatment (see [5] for a review and t54j for a 
field-theoretical approach) . 

Experimental, Monte Carlo and analytical studies of the DAFF/RFIM agree that 
this system experiences a phase transition for more than two spatial dimensions. At 
low temperatures, if the applied magnetic field h is small, the DAFF remains in an 
antiferromagnetically ordered state. If the magnetic field (or the temperature) is raised, 
however, the ferromagnetic part of the Hamiltonian weakens the antiferromagnetic 
correlations and the DAFF crosses into a paramagnetic state. This transition occurs 
across a smooth critical surface in the {j3,h,p) space (p characterizes the dilution of 
the system). 

That much is generally accepted. However, despite decades of interest on this sys- 
tem, most of the details of this phase transition remain controversial. For instance, the 
analytical work on the RFIM favors a second-order scenario. Most recent numerical and 
experimental studies of the DAFF agree |55ll56l[57ll58) . but some recent work [59II60| 
has found signatures characteristic of first-order transitions (such as metastability, or 
a very low critical exponent /? that raises the possibility of a discontinuous transition). 
These apparent inconsistencies between numerical work and analytical predictions have 
even led some researchers to consider the possibility that the RFIM and DAFF may not 
be equivalent after all. Even allowing for this possibility, the comparison of numerical 
and experimental work on the DAFF has not always yielded completely satisfactory 
results. For instance, some experiments have reported a divergent specific heat [61II62J . 
not observed numerically [56' . 

The problems on the numerical front can be understood if we notice that the DAFF 
is particularly ill-suited to numerical investigation in the canonical ensemble. On the 
one hand, the system suffers from a very severe thermally activated critical slowing 
down, where the characteristic times grow as logr ~ f , with 6 ~ 1.5. On the other 
hand, the statistical averages in the canonical ensemble are dominated by extremely 
rare events, which gives rise to strong violations of self-averaging |63j . 

Finally, some crucial physical results are hardly accessible to a canonical compu- 
tation. Perhaps the most important is the hyperscaling violations exponent 9. This is 
an additional critical exponent that one must introduce to make the theory consistent, 
and that modifies the standard hyperscaling relation, 2 — a = Dv, in the following way, 

2-a = u{D-e). (45) 

In the canonical ensemble, the exponent 9 is associated to the free-energy barrier 
between the ordered and the disordered phase 55,64 , AF oc L . Hence, it is a good 
marker to distinguish a first-order scenario {9 > D — 1) from a second-order one 
(6 < D — 1). At the same time, it is extremely difficult to compute with standard 
Monte Carlo methods. 

In [35] we performed a Tethered Monte Carlo study of this system and found 
that many of the problems with previous numerical work were solved. As one wishes 
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to ascertain that a first-order transition with the magnetic field is not present, it is 
only natural to tether the conjugate magnitude, the standard magnetization. On the 
other hand, the largest barriers are associated to the order parameter (the staggered 
magnetization, see below), which was also tethered. 

A consistent picture of the phase transition, which was shown to be continuous, 
emerged. Self-averaging was restored and thermalization greatly accelerated. We were 
able to obtain the three independent critical exponents and found a (slow) divergence of 
the specific heat, in accordance with experiments. In particular, the critical exponent 9, 
so elusive to canonical computations, was determined with some precision: 9 — 1.47(2) 
(the free-energy barriers related to 9 are immediately translated into effective-potential 
barriers, readily computed with the methods of the previous sections). 

In the remainder of this section we shall explain in some detail the application 
of the tethered formalism to this problem. We first define the model and some basic 
observables in l4.1l This is followed in 14.21 bv an examination of the DAFF with tradi- 
tional (canonical) methods, that expands on the problems discussed above. Next, we 
construct a tethered formalism adapted to the DAFF and discuss how to restore self- 
averaging to the disorder averages. In 14.51 we report the parameters of our simulations 
and in 14.61 we discuss their thermalization and possible optimizations. Finally, in 14.71 
we include a geometrical investigation of the instanton configurations, not included 
in [35]. This serves as a good example of how a tethered computation is not restricted 
to reproducing canonical results more efficiently, but can access interesting information 
that is completely hidden from a canonical study. 



4.1 Model and observables 

We consider the diluted antiferromagnet in a field (DAFF) on a square lattice with 
N — L nodes and periodic boundary conditions 

H — 2_^ ^xSx>^ySy — h2_^<^xSx. (46) 

(x.y) X 

In addition to the Ising spins [sx ~ ±1), the Hamiltonian includes occupation vari- 
ables ex- These are randomly distributed according to P(e) = p5{e — 1) + (1 — p)S{e). 
Throughout this paper we shall consider a spin concentration oi p — 0.7. This value 
is very far from the percolation threshold of p ~ 0.31 [§5], but not so high as to be 
dominated by the pure case {p = 1). The ex are quenched variables, which means 
that we should perform first the thermal average for each choice of the {ex} and only 
afterwards compute the disorder average. In the canonical ensemble, we shall denote 
the thermal average by (• ■ ■ ) and the disorder average by (■•■)■ 

For each spin configuration, we can define the standard and staggered magnetiza- 
tions 

M = Nm = ^exSx, (47) 

X 

Ms = Nnis = ^ exSxTTx, ttx = e''^^"=^ ^". (48) 
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Time (pai'allel tempering updates) 

Fig. 2 (Color online) Monte Carlo histories of the first-neighbor connectivity density u for 
different runs of the same sample in a canonical parallel tempering simulation (L = 24, T = 
1.6, h = 2.4, p = 0.7). After a variable simulation time, the runs eventually reach one of two 
local minima (labeled by ui and U2), but the simulations are not long enough for us to see 
any tunneling between the two. 



We shall distinguish between the energy E and the first-neighbors connectivity U, 

exSx, (49) 



E = Ne — 2_^ i^xSx^ySy — h2_^ ' 



tyby. 



(50) 



(^,y) 



4.2 The DAFF in canonical Monte Carlo simulations 



In this section we shall explain the limitations inherent in a canonical study and demon- 
strate how the tethered formalism can be used to gain more information about the 
system. In order to do so, we shall first carry out a comparative study of a single disor- 
der realization (sample) in both the canonical and tethered ensembles. We shall then 
explain how to perform the disorder average in the tethered formalism and provide the 
motivation for the approach taken in this study. 

We have performed several parallel-tempering [SlllU] simulations of the same disor- 
der realization for an L = 24 system|_| In Figure [2] we can see how after some time the 
systems reach one of two local minima of the effective potential, with first-neighbor 



connectivities, Eq. (|5Up 



/I) = -1.07735(8) and m^^' = -1.0973(7). The corresponding 
energy densities, Eq. ^, are e^^^ = -1.37608(2) and e^^^ = -1.3820(2). 

On the basis of these data, one would think that these two minima are relevant 
to the equilibrium of this particular sample, so that if we wait a sufficiently long 
time we will eventually see the system tunnel back and forth between the two of 
them. This apparent metastability could lead to the interpretation that the system is 



^ We used 40 evenly spaced temperatures in the range 1.6 < T < 2.575, choosing h so that 
h/T = 1.5. The dilution was p = 0.7. These parameters are taken from |60| . although in that 
reference even lower temperatures are explored. 
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1.01 1.02 1.03 1.04 1.05 1.06 1.07 1.08 1.09 1.1 1.11 
-u 

Fig. 3 (Color online) Scatter plot of the first-neighbor connectivity and the staggered magne- 
tization, as computed on the very last 0.1% of the Monte Carlo history for the simulations in 
Figure[2] Most of the simulations have reached one of two local minima with opposite staggered 
magnetization. 



experiencing a first-order phase transition. There are two problems, however. The first 
one is that we have not actually seen the tunneling in any of our simulations|f|ln other 
words, canonical parallel-tempering is not able to thermalize the system in a reasonable 
amount of time. Notice that obtaining the two minima in separate simulations is not 
enough, we need to see the tunneling in order to know their relative weight. 

The second problem is illustrated on Figure |31 In it we represent a scatter plot of 
the staggered magnetization against u for the runs in Figure [21 It is readily apparent 
that the two local minima correspond to antiferromagnetic systems with opposite sign 
of the order parameter (m-g = 0.5023(3), rris — —0.543(3)), so the observed metasta- 
bility does not correspond to the disordered-antiferromagnetic transition in which we 
are interested|j A canonical simulation cannot separate the different staggered magne- 
tization sectors, so a study in this statistical ensemble will always be contaminated by 
this spurious metastability and dominated by extremely rare events. 



4.3 The DAFF in the tethered ensemble 



Following the results of our canonical investigation, the obvious quantity to tether in 
the DAFF is the staggered magnetization mg. This would, for instance, allow us to 
study separately the different saddle points, as discussed in Section [3.21 and, through 
the Helmholtz potential, to determine their relative weight, without any need for expo- 
nentially slow tunneling. In addition, since we are interested in a transition at non-zero 
magnetic field h, in which metastability could also appear, we will tether the standard 
magnetization m. 



^ We have run 100 thermal histories taking 10^ parallel-tempering steps in each, which 
suggests that the tunneling probability is upper bounded by 10~^. 

^ Notice that an individual sample of the DAFF does not have the usual Z2 symmetry in 
antiferromagnetic systems, since the number of spins aligned with h is a. random variable even 
in a fully antiferromagnetically ordered configuration. 
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Fig. 4 (Color online) To-p: Result of a tethered simulation for the sample of Figure [2] The 
two saddle points are joined by a straight path (m(t), ms(t)), along which we measure the 
energy density {u)m(t),fn,{t)- ^"'^ t = and t = 1 we reproduce the energies of the local 
minima obtained in the canonical simulation. Bottom: Relative weight of the points along the 
simulated path, obtained from the line integral of the tethered magnetic field I I52I I (see text). 



We, therefore, have a two-component tethered magnetic field 



Vn]\f{m,rhs) — 



dm ' 9ms 

We use quadratically added demons, as in Eq. psp . Therefore, 

1/2 - 1/JV 



— (,\0)r?i,ms5 \Os)m,T?ia j (^^) 



6= 1 
fes = 1- 



m — m 

1/2 - 1/JV 
rhs — rus 



(52) 
(53) 



In general, canonical averages in the presence of an external magnetic field with 
a regular component h and a staggered component hs can be recovered through the 
Legendre transform equation: 

f f^-ffl f (j^-fil^ Q-NinN{rh,ms)-ml3h-mBfihs] ' 

although here we shall be interested in hs = 0, so we shall use the notation 

{O){h)^{O)ih,hs^0). 



{0}{h,hs) = 



(54) 



(55) 

The computation of these canonical averages following the procedure described 
in Section [2] is very complicated. It would involve the non-trivial computation of a 
two-dimensional potential from its gradient, after running a set of many simulations 
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in a two-dimensional grid in (m,ms). As we shall see in the following, the practical 
application is simpler than that. This is because for a given value of h only a very 
small region of the configuration space (m, rhs) will have a relevant weight in (|54p . as 
discussed in Section [3.21 Therefore, due to the ensemble equivalence property, we will 
be able to relate canonical and tethered averages through saddle-point equations 

(56) 
dfipf 



. drhs 



= {bs}r 



As an example we can return to the sample of Section 14.21 Using a canonical 
parallel-tempering simulation, we had identified two local minima, but were unable 
to determine their relative weights. We can perform now this study in the tethered 
formalism. 

The first step in this computation is identifying the values of {m,rhs) that cor- 
respond to the local minima using H56p . Notice that from the canonical simulation 
we know the values of m^''' and nis , so, using Eq. (|52[) . rhs — "ig -I- 1/2 and 
m'*' ~ m'*' -I- 1/(2(1 ~ j3h)). From these starting guesses the actual minima are readily 
found. The strategy is then clear: by performing the line integral of 

B = {{h)rn,ri^, " PK (&s)™,mj (57) 

along a path connecting the two saddle points we immediately obtain their potential 
difi^erence or, equivalently, the ratio between their weights. The very same procedure, 
carried out for hard-spheres, is depicted in Fig. 1151 

Figure |4] shows the result of this computation. For our connecting path we use the 
straight line {m{t),rh.s{t)) = (m^ Krhs )(1 — t) + {rh^ ,'rhs ). As we can see in the 
upper panel, the tethered averages {u)^u\ fhJt) fo'' ^ = and t — 1 closely match u^ ' 
and u^ ' , confirming ensemble equivalence. The lower panel shows the relative weight 
P(rh[t),m,s{t)) of the points along the path. In accordance with (|54p this is simply 

P(m(i),m.(i)) = ^-'N[u^{Mt),rh^{t))-l3hrh{t)\ ^ ^gg^ 

where we have chosen the zero in Q]^ so that the weight of the whole path is normalized. 
One of the two peaks is seen to have about ten times more importance than the 
other. Interestingly, we see that they are separated by a region of very low probability, 
which explains the difficulty of the canonical simulations to tunnel between the two of 
them. The large number of intermediate maxima and minima in the effective potential 
are a real effect, and not fluctuations, as we shall see more clearly in the next section. 



4.4 Self-averaging and the disorder average 

In order to perform a quantitative analysis of the DAFF, we have to simulate a large 
number of samples and perform the disorder average. The naive way of doing this for 
a system with quenched disorder would be to measure B and construct Qj^ for each 
sample, then use Eq. H54p to compute all the physically relevant {0)(h). Only then 
would we average {0){h) over the disorder. 



19 




Fig. 5 (Color online) Left: Staggered component of the tethered magnetic field, (6)^ ^^ for 
rri = 0.12 as a function of rfiB- We plot the results for several samples of an L = 24 system at 
/3 = 0.625. The curves are cubic splines interpolated from 33 simulated points. Right: Disorder- 
averaged (fes)^ ^ for the same smooth magnetization and temperature as the left panel, for 
all our system sizes. The plot shows the different behavior of the regions inside the gap, where 
the staggered magnetic field goes to zero as L increases, and outside of it, where there is a 
non-zero enveloping curve. 



This approach is, however, paved with pitfalls. First of all, computing a two- variable 
f2j^{rh, rfis) from a two-dimensional (m, rhs) grid is not an easy matter. In the previous 
section we avoided this problem by finding the local minima first and then evaluating 
^jv only along a path joining them. But this cannot be done efficiently and safely for 
a large number of samples. Even if it could be done, the free-energy landscape of each 
sample is very complicated, with many local minima, several of which could be relevant 
to the problem, so a very high resolution would be needed on the simulation grid. 

Finally, even reliably and efficiently computing the canonical averages {0){h) would 
not be the end of our problems. Indeed, it is a well known fact that random systems 
often suffer from violations of self- averaging [66ll67lf55| . This phenomenon has recently 
been studied in detail for the RFIM [63], where the violations of self- averaging are 
shown to be especially severe. In this section we address all these problems and demon- 
strate the computational strategy followed in our study of the DAFF. 

The first step is ascertaining whether the tethered averages themselves are self- 
averaging or not. We already know that, since we are going to be working with a large 
regular external field h, the relevant region for the regular magnetization m (and hence 
m) is going to be very narrow. Therefore, we are going to explore the whole range of 
rhs for a fixed value of m = 0.12. 

We have plotted the staggered tethered magnetic field {bs)m=o.l2,ms foi' 20 samples 
of an L = 24 system at /3 = 0.625 in Figure[5] — left. The different curves have a variable 
number of zeros, but all of them have at least three: one in the central region and two 
roughly symmetrical ones for large staggered magnetization. The positions of the two 
outermost zeros clearly separate two differently behaved regions. Inside of the gap the 
sample-to-sample fluctuations are chaotic, while outside of it the sheaf of curves even 
seems to have an envelope. This impression is confirmed in the right panel of Figure O 
where we show the sample-averaged tethered magnetic field for several system sizes. 
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Fig. 6 (Color online) Left: Sample-to-saniple fluctuations in the staggered component of the 
tethered magnetic field, Eq. II59II . against system size. We show two sets of points, both at 
/3 = 0.625 and m = 0.12, but at different values of rria- The red curve, corresponding to 
ms = 0.7 is right inside the gap defined by the outermost zeros in Figure [S] while the blue 
curve (ma = 1.1) is outside. Both are shown to decay with a power law. Right: Disorder average 



(6s>, 



for the systems of the left panel. Inside the gap, the field goes to zero with L; outside 



it has a finite limit. 



In order to quantify this observation we can study the fluctuations of the disorder- 
averaged {6s)a,a , 



'(m,ms) = I (is) 



{bs}r 



(59) 



This quantity is plotted on the left panel of Figure [^ As we can see, it goes to zero as 
a power in L, so we also plot fits to 



CT (m = 0.12, 771s) = A(rhs)L 



c^rhs) 



(60) 



This would seem like a very good sign, because it could be indicative of self-averaging 
behavior. However, it is not the whole story. If we recall the right panel of Figure [SI 
we see that inside the gap the tethered magnetic field itself, not only its fluctuations, 
goes to zero as L increases. In fact, as shown in Figure [6] for rhs = 0.7, a ~ L~ , 
while (fes),^ jji ~ L^ ■ This means that the relative fluctuations o-/ (bs)^ ^ do not 
decrease with increasing L. For the point outside the gap, however, the disorder average 
of the tethered magnetic field reaches a plateau. Furthermore, the fluctuations decay 
with c = 2.94(7), which is compatible with the value c = D that one would expect in 
a self- averaging system. 

In physical terms, this analysis means that the local minima for a small, but non- 
zero, value of the applied staggered magnetic fleld hs would be self-averaging. We can 
now recall the well-known recipe for dealing with spontaneous symmetry breaking: 
consider a small applied field and take the thermodynamical limit before making the 
field go to zero. Translated to the DAFF, this means that we should solve the saddle- 
point equations (|56p on average, rather than sample by sample, and then take the 
/is — >■ limit on the results. 






(61) 



drh: 
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Fig. 7 (Color online) Disorder-averaged probability density function of rhs conditioned to 
m = 0.12 for the system of Figure [5] 



In other words, we are considering the disorder average of a thermodynamical potential, 
f2j\j, different from the free energy. This approach was first introduced in [B5], in a 
microcanonical context (the averaged potential was in that case the entropy). 

Notice that for the regular component of the tethered magnetic field we are always 
going to work with a finite value of the external field h. Therefore, the integral in (|54p 
is going to be completely dominated by an extremely narrow m range. Therefore, 
we can consider different values of rfi separately and relate them to the particular 
h that generates a local minimum at that magnetization. This is best accomplished 
by integrating the disorder-averaged tethered magnetic field along the path (m = 
mo,rhs(t)). In this way we obtain the probability distribution of rhs, conditioned to 
m — iriQ, which we will denote by p(ms|mo) (Figure [T}. This probability density 
function can be used to average over rhs for fixed m. 



(O), 



/ dms p 



where 



(ms|m)(0). 



dnN{ms\iri) 
drhs 



I 



drhs (O),, - e" 



-N f2 N (^s 1^) 



{bs)r 



(62) 



(63) 



With this process we have integrated out the dependence on rhs. We now have a 
series of smooth functions (O)^, together with a smooth one-to-one function h{rh) — 

(^)m/P- Remembering our discussion of ensemble equivalence in Section [3.21 we shall 
make the approximation 

(Oy(/i)^M™(h)- (64) 

At any rate, we remark that the r.h.s. is better suited for reproducing the physics of 
experimental samples. 

Finally, let us mention that some of the physical observables are correlated with 
the distribution of the ex . This can be exploited to reduce the statistical errors of their 
sample averages, using the technique of control variates [701140] . 
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Table 2 Parameters of our simulations. For each of the A^aamples disorder realizations for 
each L, we run A^^ X Af^^ tethered simulations with temperature parallel tempering. The 
Nt participating temperatures are evenly spaced in the interval [1.6,2.575]. For each size we 
report the minimum number of Monte Carlo steps (Metropolis sweep + parallel tempering). 
After the application of our thermalization criteria, some of the simulations for L > 24 needed 
to be extended, leading to a higher average number of Monte Carlo steps {N^q). 



L 


-^^samplcs 


Nt 


Nrh 


Nrh, 


ATmin 

"mo 


IV av 

"mo 


8 


1000 


20 


5 


31 


7.7 X 10^ 


7.7 X lO'' 


12 


1000 


20 


5 


35 


7.7 X 10^ 


7.7 X 10^ 


16 


1000 


20 


5 


35 


7.7 X 10^ 


7.7 X 10^ 


24 


1000 


40 


5 


33 


7.7 X 105 


9.3 X 105 


32 


700 


40 


4 


25 


1.5 X lO'' 


5.5 X 10^ 



4.5 Our simulations 

We can infer several useful conclusions from the analysis of the previous section 

— The disorder average should be performed on the tethered observables, before com- 
puting the effective potential. 

— It is best to analyze several values of rh separately, since the average over rhs for each 
fixed rh can be unambiguously related to the canonical average via h{rh) = {&),f,//3. 
In this way, we can study the phase transition that arises by varying the applied 
magnetic field at fixed /3. 

— For fixed m the conditioned probability p(rhs\rh) has two narrow, symmetric peaks, 
separated by a region with extremely low probability. 

Therefore, we have carried out the following steps 

1. Select an appropriate grid of m values. This should be wide enough to include 
the critical point for the simulation temperature, and fine enough to detect the 
fluctuations of (O)^. These turn out to be very smooth functions of rh, so a few 
values of this parameter suffice i39,. 

2. For each value of rh, select an appropriate grid of rhs- We start with evenly spaced 
points and after a first analysis add more values of rhs in the neighborhood of the 
minima, as this is the more delicate and relevant region. 

3. The simulations for each (rh, rhs) are carried out with the Metropolis update scheme 
of section ETTI In addition, we use parallel tempering (attempting to exchange con- 
figurations at neighboring temperatures). This is not needed in order to thermalize 
the system for L < 32, but it is convenient since we also study the temperature de- 
pendence of some observables. Furthermore, the use of parallel-tempering provides 
a reliable thermalization check (see the Appendix, particularly section |A.1|) . 

From Figure [7] it seems that, in addition, it would pay to restrict the simulations to 
a narrow range of rhs around the peaks, since the tethered values for rhs away from 
them cannot possibly contribute to the average. This is true but for one exception, the 
computation of the hyperscaling violation exponent 9, which requires that the whole 
range be explored (see [39])l I 

The parameters of our simulations are presented on Table [5] The table lists the 
number A*'^ of values in our rh grid and the number of points in the rhs grid for each, 
so the total number of tethered simulations for each sample is A^^ x N,fi^. 



* In Section l4.7l we show an additional study that requires simulations far from the peaks. 
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The number of Monte Carlo steps in each tethered simulation is adapted to the 
autocorrelation time (see next section), the table lists the minimum length A'JJ™ and 
the average length N^q for each lattice size. 



4.6 Thermalization and metastability in the tethered simulations 

We have assessed the thermalization of our simulations through the autocorrelation 
times of the system, computed from the analysis of the parallel tempering dynamics 
(see Section lA.ip . We require a simulation time longer than lOOrintU As a fail-safe 
mechanism, we compute the time t^^^^ that each configuration spends in the upper 
half of the temperature range. If any of these is less than a third of the median ihot. 
we double the simulation time, considering that the simulation is too short to allow 
a good determination of ri„f This process is only followed for L > 24. For smaller 
sizes we have simply made the minimum simulation time large enough to thermalize 
all samples. 

The distribution of correlation times for our different samples turns out to be 
dependent on the value of (m, rris). Considering first the variation of the average ri„t 
with rfis at fixed m, we see that the minimum and its adjoining region are much 
easier to thermalize (Figure (8)1 . This region coincides with the only points that have 
a non-negligible probability density (Figure (Tjl, i.e., the only points that contribute to 
the computation of the {0),j,. This fact suggests a possible optimization, that we will 
discuss in Section [4.6.1l 

A second interesting result comes from studying the evolution of the T[nt with m. 
Figure [9] represents the histogram of autocorrelation times for rhs = 0.8 (in the 'hard' 
region) for two values of m. The distribution for m = 0.12 has a much heavier tail. It 
is shown in 39 that this is due to the onset of a phase transition. 

The difficulty in thermalizing some samples stems from the coexistence of several 
metastable states, even for fixed (m, rhs). In Figure [TOl we represent the time evolution 
of the energy u for several temperatures of the same sample (L — 32, m = 0.12, 
rhs = 0.8). As we can see, for a narrow temperature range several metastable states 
compete. This has a very damaging effect on the parallel tempering dynamics, that 
get stuck whenever a configuration that is metastable for one temperature is very 
improbable in the next (see Figure [11}. 

For L = 32, some point^_j presented a metastability so severe that enforcing a sim- 
ulation time longer than 100ri„t would require a simulation of more than 10 parallel- 
tempering updates (one thousand times longer than our minimum simulation of ~ 10 
steps). Thermalizing these points (which constitute about 0.1% of the total) would 
have thus required some 10 extra CPU hours, with a wall clock of many months. 
We considered this to be disproportionate to their physical relevance (they are all re- 
stricted to a region far from the peaks where the probability density is < 10~ , see 
Figure [T]). Therefore, we have stopped these simulations at about lOrint. This is still 
a more demanding thermalization criterion than is usual for disordered systems and 



^ For most simulations Ti^t ~ fcxp, so this value is mucli larger than what is required to 
achieve thermalization. This ample choice of minimum simulation time protects us from the 
few cases where Tcxp is noticeably larger than Tint- 

^^ By 'point' we mean any of the A^samplcs x ^m x ^rhs individual tethered simulations for 
each L. 
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Fig. 8 (Color online) Histograms of thermalization times for our L = 32 simulations, at two 
values of rhs for m = 0.12. Notice the logarithmic scale in the horizontal axis, which is in 
units of 300 parallel-tempering steps. The minimum of the effective potential (the peak in the 
probability distribution) at this fh = 0.12 corresponds to riis = 1.078. Notice that we cannot 
measure times shorter than our measuring frequency of 300 parallel-tempering steps, so the 
first bin should be taken as encompassing all shorter autocorrelation times. Only the samples 
with T > 2*' have to be extended from the minimum simulation time. 



does not introduce any measurable bias in the physically relevant disorder-averaged 
observables. This can be checked in several ways: 

— First of all, as we have already discussed, these points are restricted to a region in 
TTis with probability density of at most 10~ . Therefore, even if there were a bias 
it would not have any effect in the computation of canonical averages. 

— For the only affected physical observable, the free-energy barriers used to compute 
9 (see [31]), we can compare the result for L = 32 with the extrapolation from 
smaller sizes, and it is compatible. 

— Even at the most difficult values of {m,rhs) the log2-binning plot (the only ther- 
malization test typically used in disordered-systems simulations) presents many 
logarithmic bins of stability (Figure [T^] — left). Even if we subtract the result of 
the last bin from the others (in equilibrium, this substration should yield zero), 
taking into account statistical correlations [71], several bins of stability remain 
(Figure [T2] — right). This is a very strict test, and one that even goes beyond phys- 
ical relevance (because it reduces the errors dramatically from those given in the 
final results). 

4-6.1 Optimizing Tethered Monte Carlo simulations 



We can highlight two interesting facts from the previous discussion 

— Only a very narrow region around the minima of the effective potential has any 
significant weight for reconstructing canonical averages (Figure [7]). 

— It is much harder to equilibrate the region far from the minima (Figure [Sjl. 

Combining these two observations, it turns out that, if our only interest is reconstruct- 
ing canonical averages, we can achieve a qualitative improvement in simulation time 
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Fig. 9 (Color online) Same as Fig. [S] but now we consider several values of m for rhs = 0.8 
(in the hard therinalization region far from the peak). The points for rh = 0.12, closer to the 
critical point, exhibit a much heavier long-times tail (mind the vertical logarithmic scale). 
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Fig. 10 (Color online) Time evolution (in number of parallel tempering steps) for the energy 
density u for several temperatures of a single L = 32 sample at rh. = 0.12, rhs = 0.7. For most 
temperatures the equilibrium value is quickly reached, but for a very narrow temperature range 
there are several competing metastable states (bottom-left panel for T = 1.901). 



by simulating only a narrow range around the peak in P(ms|m) for each value of the 
smooth magnetization m. In fact, all the results discussed in |39j could have been 
computed in this simplified way, except for the hyperscaling violations exponent 9. 

We have demonstrated this optimization by simulating 400 samples of an L = 48 
system for rh = 0.12. We use only A'^^ — 6 (but we have to increase the number of 
temperatures in the parallel-tempering to keep the exchange acceptance high). Table [3] 
shows the value of mS°^ as a function of L. 
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Fig. 11 (Color online) Temperature random walk for one of the 40 replicated systems for 
the parallel tempering simulation of the sample in Figure 1101 The flow is blocked at the same 
temperature that had several metastablc states (marked with a horizontal line) , see bottom-left 
panel in Fig. 1101 
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Fig. 12 (Color online) Left: Time evolution of (fes)^ „, for rh = 0.12, rrig = 0.5, in a logarith- 
mic scale (the first bin is the average over the second half of the simulation, the second bin the 
average over the second quarter and so on). This is in the middle of the (m,ms) region where 
thermalization is hardest. Right: Same figure, but now subtracting the value of the first bin 
from each of the following ones. Since the data are correlated, this subtraction causes a very 
significant error reduction, but the difference is still compatible with zero for several bins. 



Table 3 Value of the peak position for our different system sizes. The result for L = 48 is 
of comparable accuracy, despite being computed from only 6 tethered simulations around the 
minimum of the effective potential. Notice that for L < 32 we can average over the positive 
and negative peaks, so the number of samples for these systems is effectively double the value 
shown in Table (2] 



K 



samples 



Nr^„ 



■ peak 



1/2 



0.58585(87) 
0.58239(54) 
0.58154(36) 
0.57839(24) 
0.57672(20) 
0.57491(33) 



8 
12 
16 
24 
32 
48 



1000 X 2 
1000 X 2 
1000 X 2 
1000 X 2 
700 X 2 
400 
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Fig. 13 (Color online) Equilibrium configuration for an L = 24 system at /? = 0.625, rh = 0.12, 
ms = 0.5. 



4.7 Geometrical study of the critical configurations 



The picture painted in [35] is that of a second-order transition, but one with quite 
extreme behavior. Among its peculiarities we can cite an extremely small value of j3 
and large free-energy barriers, features both that are reminiscent of first-order behavior. 
The free-energy barriers in the DAFF, however, diverge too slowly to be associated to 
the surface tension of well defined, stable coexisting phases. But this only begs the 
question of what kind of configurations can give rise to such a behavior. 

In this section we shall study the geometrical properties of the minimal-cost spin 
configurations joining the two ordered phases at the critical point. To this end, we 
consider simulations aX j3 = 0.625, m = 0.12 « wicritical ^-nd ms = 0.5. Recalling that 
ms — 7718-1-1/2, this last condition expresses the fact that we are studying configurations 
with no global staggered magnetization. This is a good example of an 'inherently 
tethered' study, that examines information hidden from a canonical treatment. 

Figure [13] shows an example of such a configuration for an L — 24 system. In order 
to make the different phases clearer, we are not representing the spin field Sx, but the 
staggered field SxTTx. As is readily seen, even if the global magnetization is ms ~ 0, 
the system is divided into two phases with opposite (staggered) spin. In geometrical 
terms, most of the occupied nodes of the system belong to one of two large clusters 
with opposite sign, with only a few scattered smaller clusters. 

At a first glance, this picture may seem consistent with a first-order scenario, were 
the system is divided into two strips whose surface tension gives rise to the free-energy 
barriers (see [31) for an example of these geometrical transitions in a first-order setting). 
In order to test this possibility, we can study the evolution of the interface mass with 
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Table 4 Masses of the interfaces for our equilibrium spin configurations at /3 = 0.625, m = 
0.12 Si iric, rhs = 0.5. This mass grows as Afintorfaco <^ L'^- From a fit, c = 2.2402(24), which is 
incompatible with our estimate of 8 = 1.469(20), confirming that the free-energy barriers are 
not associated to surface tensions. 



L 


-^'' samples 


-'^interface 


12 


3000 


160.29(25) 


16 


3000 


304.62(41) 


24 


3000 


755.74(94) 


32 


700 


1446.1(39) 



the system size and compare it with the explicit computation of free-energy barriers 
done in [39) . 

Given a configuration, we first trace all the geometric antiferromagnetic clusters. 
We then identify the largest and second largest ones. Finally, we say that an occupied 
node belongs to the 'interface' if it belongs to the largest cluster and has at least one 
first neighbor belonging to the second largest one. We have computed in this way the 
interface A^intortacc for our 700 L = 32 samples and for 3000 samples for all our smaller 
systems (we have run additional simulations just at this point). Table [4] shows the 
result of this computation. A fit to 

A^intcrfacc = ^L" , (65) 

for L > 12 gives c = 2.2402(24) with x^/d.o.f. = 3.62/2. One may think that such 
a large exponent could be indicative of a surface tension. However, the study of |39| 
showed that free-energy barriers in the DAFF grow as L , with 6 = 1.469(20) 7^ c. 
Indeed, for the archetypical second-order model, the D-dimensional pure Ising ferro- 
magnet, the percolating clusters are space filhng, so c = D. 

A second interesting feature of the configuration pictured on Figure [TS] is that 
there is a path connecting the spins in the green strip across the red one (there is 
only one green strip, since we are considering periodic boundary conditions). If we 
were to study a complete tomography of this configuration, we would find several of 
these paths (which, of course, need not be contained in a plane). In other words, the 
phases are porous. This is in clear contrast to a first-order scenario in which the phases 
are essentially impenetrable walls. Now, this could be a peculiarity of the particular 
selected configuration. In order to make the analysis quantitative, we shall examine 
all of our samples and determine the strip-crossing probability P. This is defined as 
the probability of finding a complete path with constant staggered spin across the 
strip with opposite staggered magnetization and we can compute it with the following 
algorithm: 

1. For each configuration, compute the Fourier transform of the staggered spin field 
at the smallest nonzero momentum in each of the three axes ((px, (py, (pz)- 

2. If the system is in a strip configuration, one of the (p will be much larger than the 
other two. Assume this is (px, so the strips are perpendicular to the OX axis. 

3. Measure the staggered magnetization A/^ on each of the planes with constant x 
and identify the plane x = a^max with largest |M^|. This plane will be at the core 
of one the strips. 

4. Trace all the clusters that contain at least one spin on the x — a^max plane, but 
severing the links between planes x = Smax and x — a;max — 1. 




Fig. 14 (Color online) Probability of completing a path with constant staggered spin across 
the strip with opposite staggered magnetization, see Fig. 1131 for our spin configurations at 
/3 = 0.625, rhs = 0.5. 



5. If any of the clusters reaches the plane x = Xmax — 1 there is at least a path through 
the strip with opposite magnetization (the previous step has forced us to go the 
long way around, so we know we have crossed the strip). 

We have plotted the strip-crossing probability P as a function of rh in Figure 1141 
Notice that 1 — P behaves as an order parameter. If we keep increasing rh, so that 
we enter the ordered phase, the phases eventually become proper impenetrable strips, 
hence P = for large enough systems. On the other hand, for low m, in the disordered 
phase, the strips become increasingly porous, so that P = 1 in the limit of large 
systems. In fact, the inversion of the finite-size evolution at m ~ 0.12 signals the onset 
of the phase transition. 



5 Hard Spheres Crystallization 

In this Section we consider a completely different problem, namely hard-spheres crystal- 
lization, as an illustration of the flexibility of the Tethered Monte Carlo. In particular, 
this problem has neither quenched disorder, nor a supporting lattice. Furthermore, the 
phase transition is of the first-order, rather than continuous. 

Actually, crystallization (or melting) is probably the most familiar example of a first 
order transition. Most fluids undergo a transition to a crystalline solid upon cooling 
or compression. One might thus be shocked to learn that a bona-flde Monte Carlo 
simulation of crystallization, conforming to the standards taught in good textbooks [45] . 
is plainly impossible nowadays. 

The problem is that there are simply too many local minima of the effective po- 
tential where the simulation can get trapped. Some degree of cheating (you could also 
say artistry) is needed to guide the simulation. 

That the problem is a fairly subtle one is illustrated by the fact that hard-spheres 
crystallize in three dimensions [721173) . even if the fluid and the face-centered cubic 
(FCC) crystals have exactly the same energy. In fact, hard spheres have become the 
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standard model in the field, where all new ideas must be tested. Besides their theoretical 
interest, hard spheres provide as well an important model for colloidal suspensions |74l 

Eg. 

In fact, workable numerical methods [76j put by hand the crystal in the simulation. 
One may try to achieve equilibrium between the crystal and the fluid, as in the phase 
switch Monte Carlo ^77,, but this is feasible only for small systems (up to A'^ = 500 hard 
spheres [7S]). One may also compute separately the fluid and solid free energies. For the 
fluid, one resorts to thermodynamic integration, while several methods are available 
for the crystal (Wigner-Seitz [79], Einstein crystal [SOllST] . Einstein molecule [82]). The 
transition point is determined by imposing to both phases the conditions of equal 
pressure, temperature and chemical potential. An alternative is direct coexistence [831 
184] . a dynamic, non-equilibrium method that simulates rather large systems with great 
accuracy [85] . 

Our scope here is to illustrate how Tethered Monte Carlo can be used in this con- 
text. Indeed, it has been recently shown that Tethered Monte Carlo allows equilibrating 
up to A'^ = 2916 hard spheres at their phase coexistence pressure I38J . As explained in 
the Introduction, constrained Monte Carlo studies of crystallization kinetics have been 
carried out before [3]|33). in an umbrella sampling framework. 

The remaining part of this section is organized as follows. In Sect. l5.l1 we recall the 
hard sphere model. Our order parameters are defined in Sect. 15.21 The specific imple- 
mentation of the tethered formalism is in Sect. 15.31 Details about the computation of 
the phase-coexistence pressure 38, are provided in Sects. [5^ and [531 The performance 
of the Tethered Monte Carlo algorithm is assessed in 15.61 Finally, we comment on the 
computation of the interfacial free energy in Sect. 15.71 



5.1 The hard spheres model 

We consider a collection of A'' hard spheres, of diameter o. They are contained in 
a cubic simulation box, with periodic boundary conditions. The system is held at 
constant pressure p (hence the simulation box may change its volume, but remaining 
always cubic). 

Let us introduce the shorthand R for the set of particle positions, {ri}^^^. The 
constraint of no overlapping spheres is expressed with function H[R), which vanishes 
if any pair of spheres overlaps {H{R) = 1 otherwise) I I 

The Gibbs free-energy density, g{p,T), is obtained from the partition function 



Y^,^ = e-^'^^(-^) = -^ l^ dl/e-^''^ I dRHiR) . 



where A is the de Broglie thermal wavelength, while /? — l/{kQT). 



(66) 



^^ For a standard model fluid, one would replace H{R) by a Boltzmann factor 
exp[-U(R)/ikjiT)] . 



31 
5.2 The two order parameters 

The standard order parameter for modern crystaUization studies is Qe |34II86] , which 
is the I = 6 instance of 



(67) 




where 






Qim = ^^' ^'"'^ ^ , qUi) = V YUr,,) ■ (68) 

\ ^ -^ " AT. t^\ ^ * 



e:liA^6« 



j=i 



In the above expression, YimiTij) are the spherical harmonics, while f^ is the unit 
vector in the direction joining particles i and j, rij = r, — Vj. Ni,{i) represents the 
number of neighbors of the ith particle. |_J It is important to note that, setting aside 
the periodic boundary conditions for our simulation box, Qg is rotationally invariant. 
Qq is of order 1/\/N in a fluid phase, whereas Qg ~ 0.574 or 0.510 in perfect FCC or 
BCC crystals, respectively. However, defective crystals have lesser values (Qg ~ 0.4 is 
fairly common). 

Now, if one tries to perform a tethered computation, choosing Qq as order param- 
eter, it is soon discovered that it is almost impossible to form an FCC crystal if the 
starting particle configuration is disordered. Instead, the simulation gets stuck in heli- 
coidal crystals, with a similar value of Qg, which are allowed by the periodic boundary 
conditions. The crystalline planes for these helicoidal structures are misaligned with 
the simulation box. Upon reflection, one realizes that the problem lies in the rotational 
invariance of Qg. When a crystalline grain starts to nucleate from the fluid, which is 
encouraged by the tethered algorithm for large Qg, we face the problem that Qg con- 
veys no information about the relative orientation of the growing grain with respect 
to the simulation box. At some point, the misaligned crystalline grain is large enough 
to hit itself through the periodic boundary conditions, and it needs to contort to get 
some matching for the crystalline planes. 

The way out is in an order parameter with only cubic symmetry. Such a parameter 
was recently proposed |87) : 



2288Lj=iLj=i Ca(r-,j) 64 



C-^-^ —'-rJ-,,: ' ~^, (69) 

where 



79 El^m^ 79 



Ca(f-) = x^y\l - /) + x^z\l - /) + y^z\l - x^) . (70) 

The expectation value for C in the different phases is the following: 0.0 in the fluid, 1.0 
in the ideal FCC crystal, perfectly aligned with the simulation box, and —0.26 in the 
perfectly aligned ideal BCC. The difference with the quoted value in Ref [87] for the 
perfect BCC crystal is due to our smaller threshold for neighboring particles. Again, 
as it happens for Qg, we must expect lower \C\ values for defective structures. 



^^ Two particles i and j are considered neighbors iff Tij < 1.5 cr. This choice ensures that 
we enclose only the first-neighbors shell in the FCC structure, for all the densities of interest 



here 111] . 
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In fact, we find that tethered simulations with a single order parameter (C) and 
large C, equilibrate easily. One forms nice crystals even if the starting particle configu- 
ration is disordered, and the Monte Carlo expectation values turn out to be independent 
of the starting configuration, as they should. 

However, not all is well. At some intermediate values of C the simulation suffers 
from metastabilities. Heterogeneous states with a slab of FCC crystal in a fluid ma- 
trix are degenerate with helicoidal crystals, that fill most of the simulation box, but 
have a small value of C because of their misalignment. Fortunately, these two types 
of competing configurations are clearly differentiated by Qg, which is large for the he- 
licoidal crystal, but small for the crystalline slab surrounded by fiuid. Therefore, the 
combination of the two order parameters labels unambiguously the intermediate states 
between the fiuid and the FCC crystal. 



5.3 Tethered formalism for a hard sphere system 

Since we wish to constraint simultaneously the two crystal order parameters, Qg and 
C, we use the formalism in Sect. 13.11 We choose to couple the demons linearly, which 
frees us from the constraints Qq > Qq{R) and C > C{R) imposed by quadratic 
demons. Note that ascertaining thermalization is an issue in crystallization studies. It 
is very important to compare the outcome of simulations with widely differing starting 
configurations. In this respect, the constraints Qg > Q(){R) and C > C{R) are a major 
problem, as they prevent us from using the ideal FCC crystal as starting configuration. 

As for the tunable parameter a, see Eq. (|37|) . we choose a = 200. In this way, the 
convolution in the probability distribution functions induced by the demons does not 
result into a gross distortion! I 

Our tethered statistical weight is 



u^^iR, V;p, Qg, C) = ^HiR) ,-f^r>V-N. [{Q,-Q,(R)f+[c-C(R)f] /2 
As follows from Eq. H38I) . the gradient of the Helmholtz effective potential is 



(71) 



V 
with 



«jv(06,C7,p) - (^ ^^- , ^ J - [y>Q.)Q„c,p^\^c)Q^^c,p) (^2) 



hq, =a{Qe~ Qg) , be = a {C ~ C) . (73) 



The relationship between the effective potential, 0]^{Qq, C ,p), and the Gibbs free- 
energy density is straightforward: 

FjVpT = e-^3"(P'^) = /"dQg dC e-^^^CS^'^'rf . (74) 

Furthermore, a saddle-point approximation, see Sect. 13.21 shows that 

gN{p,T) = nN{Ql,C\p) + 0{l/N), (75) 



^^ Recall that we are convolving in Eq. I I22I1 with a Gaussian of width l/\/aN. The number 
of spheres that can be equilibrated (A'^ ~ 3000) is rather modest as compared to spin models, 
which suggests enlarging a. 
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Fig. 15 (Color online) Vi7iv=256i as computed from Eq. II73II for A'^ = 256 hard spheres, 
along the straight path that joins the fluid and the FCC minima of the effective potential. 
The simulation pressure is the phase-coexistence one. To improve visibility, we have divided 
Vi7jv=256 t>y a factor of 10. Inset: map of the gradient field Vi7jv=256i including the points 
corresponding to the fluid, FCC and BCC potential minima. For the sake of visibility we have 
divided the gradient by a factor a = 200 (mind the different normalization as compared with 
main panel). 



where {Qg, C* ,p) corresponds to the p-dependent absolute minimum of Qj^{Qq, C,p), 
regarded as a function of Qq and C However, it is important to keep in mind that near 
the phase transition two local minima (namely the fluid and the FCC crystal) become 
exactly degenerate. 



5.4 The coexistence pressure: computing differences in the effective potential 



The coexistence pressure, pco, follows from the difference in effective potential between 
the pressure-dependent coordinates of the coexisting pure phases: 



AFCC 



fluid/ 



^fluid / 



An^ip) = nN{Q'e''''{p),C'^'^{p),p) - r2^(Qr^(p),C""'^(p),p) 



(76) 



^N 



The scope of the game is finding the coexistence pressure, Pco, such that Z\J7jv(pco) = 0. 
Indeed, the saddle-point condition (jTSp . tells us that, at p^o, the chemical potential for 
the two phases coincides. 

Now, the computation of Zii?7v(p) is strictly analogous to that of section [4^ We 
obtain the difference in the effective potential as a line integral of the conservative field 
B. The most convenient path is a straight line, see Fig. 1151 As in Sect. 14.31 we divide 
the path in a mesh, perform independent simulations at each point of the grid, and 
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compute AQ]^{p) from a numerical integration of the gradient, Eq. (|72[l . projected 
over the path. 

Let us stress only the main difference with the computation in Sect. 14.31 



— Here, we wish to compute Zifi^vCp) ^-s a function of pressure. Instead, in section [4^ 
we were restricted to a single value of the control parameter (temperature in that 
case). 

— The disorder-averaging procedure introduced in Sect. [43] dispensed us from a care- 
ful search of local minima. Instead, locating with very high accuracy the extremal 
points (QfC°(p),C'^CC(p)) and (Qe"''^^), C'""''^(p)) is a must here, see Sect. [531 

— The choice of a straight line as integration path, see Fig. 1151 is not only dictated 
by simplicity. It is also a matter of computational convenience. Thermalization is 
easier to achieve over it. Furthermore, as the gradient map in the inset of Fig. [15] 
shows, the gradient takes its smallest absolute value over this path. 

In order to compute Pcoi let us neglect the pressure dependence of the end points 
for the integration path in Fig. 1151 One may easily correct for end-points displacements, 
as explained in Sect. 15.51 which induces a correction in Pco negligible with respect to 
our statistical errors. 

Under the above simplifying assumption, we only need the pressure dependence of 
the gradient field, B, over the straight path in Fig. 1151 Using histogram reweighting [881 
189] , we extrapolate our numerical results at pressure p to a neighboring p -|- 5p: 

/£>\ _ J 'Q6,C,P fy.^-. 

^ 'Q6,C,p+8p- u-SpV\ . . • ^"> 

A standard argument tells us that the maximum safe extrapolation, 5p, is determined 
by the probability distribution function (pdf) for the specific-volume, v = V/N, 

c- inaxiinum t ^\^) at-P/ 2\ / \2 1 /na\ 

Hence, it is crucial that the pdf for v be unimodal (i.e. single-peaked), and with an A'^- 
independent Xp, for all point along the integration path. In other words, it is important 
that the integration path be free of metastabilities. Since this condition holds [41], it 
is straightforward to compute Ani\[{p + 5p) from simulations at p. 

Following these guidelines, p^o was computed in Ref. [5S] for 108 < N < 2916. The 
large N extrapolation was 

pgS = 11.5727(10) M, . 

The best previous equilibrium estimate seems to be the rather crude p^ = 11.50(9) [77], 
obtained using phase-switch Monte Carlo. In fact, the only previous method accurate 
enough to provide a meaningful comparison is the non-equilibrium direct-coexistence: 
Pco ~ 11.576(6) [85]- Note, however, that in order to achieve such a small error (but 
still six times larger than the tethered error), systems with up to A = 1.6 x 10 particles 
were simulated l85l. 
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5.5 Calculation of the extremal points and corrections 

We need to locate the two extremal points in the straight path in Fig. 1151 which 
correspond to the fluid or to the FCC crystal. The two points are local minima of fij^, 
regarded as a function of Qg and C but at fixed pressure. Our procedure has been as 
follows. 

We first obtain a crude estimate from standard simulations in the NpT ensemble 
(without any constrain in the crystal parameters). Note that the autocorrelation time 
for such simulations is unknown, but larger than any simulation performed to date. 
Hence, these standard simulations get stuck at the local minimum of Qj^ which is most 
similar to their starting configuration. Starting the simulation either from an ideal gas, 
or from a perfect FCC crystal, we approach the pure-phases we are interested in. The 
Monte Carlo average of Qq{F{) and C{R) provides our first guess. 

To refine the search of either of the two local minima (QqjC*), we note that, up 
to terms of third order in Qq — Qq or C — C* , 

nN{Q6,C) = n*N + ^{Q6-Qlf+AQc{Q6~Q6){C^C*) + ^iC-C*f. (79) 

The shorthand i7^ stands for n]\j{QQ,C*). Incidentally, Eq. ()79p tells us that the 
computation in Sect. 15.41 is intrinsically stable. An error of order e in the location of 
{Qq, C*) will result in an error of order e in the coexistence pressure. 

Yet, the tethered computation does not give us access to I^^Vi but to its gradient: 

Vr2^(Q6,C) = {AQQiQe-Ql)+AQc{C^C*),Acc{C-C*)+AQc{Q6~Ql)) ■ (80) 



Eq. (|80p holds up to corrections quadratic in Qq — Qq or C — C . We thus compute 
the expectation value of the field B, in a grid of nine points {Qq, C) that surround our 
first guess for {Qq, C*), and fit the results to Eq. (|80|l . We iterate this procedure until 
an accuracy ~ 10~ in both coordinates {Qq,C*) is reached [41| . 

Actually, Eq. (jTTp . shows how one extrapolates the expectation values for the gra- 
dient field from the simulated pressure, p to a nearby p + 5p. The corresponding fit to 
Eq. (jSOp provides the new coordinates [Qq{p + Sp), C* {p + 5p)\. 

At this point, one could worry because the integration path in Fig. [15] is no longer 
appropriate at pressure p -\- Sp. In fact, the extremal points in the integration path are 
pressure-dependent. However, some reflection shows that this is not a real problem. In 
fact, 

Z\J7jv(p + Sp) = A^^'^{p,p + 5p) + AP^'^{p,p + 5p) - Zi*^"''^(p,p + 5p) . (81) 

The different pieces are 

Zi''^^ (p, P + <5p) = X2^ (qI^C (p + 5p) , C^CC {p + 5p)-p + Sp) - (82) 

-nN{Qr'^{p),C^'^'^{p);p + Sp), 

the correction due to the shift of order Sp in the coordinates of the FCC minimum, 

Z\P^*^ (p, p + 5p) = r?^ (QfCC (p) , (7^^^ {p);p + 5p) - (83) 

- nN{QQ'^"^{p),C ""^(p); p + Sp) , 
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Fig. 16 (Color online) Normalized time autocorrelation function, Eq. 1871 for the specific 
volume (left) and the crystal order parameter (right), as computed for a system of N hard 
spheres, in the fluid minimum of the effective potential (labeled 5 = 0). Time is measured 
in units of EMCS (see text). Mind the different time scale for the left and right panels. Each 
value of N was simulated very close to (but not precisely at) its phase-coexistence pressure 



the line-integral sketched in Fig. \T^ as computed at pressure p + 5p, and 



.fluid/ 



(p,p + 5p) = i7jv(Q6"' (p-l-5p),C"' {p + 5p);p + 5p) 
- ^n{QT^'^{p),C ""^(p); p + 5p) , 



(84) 



the correction due to the shift in the coordinates of the fluid minimum. 

Now, one expects that the pressure-induced changes in the minima coordinates as 
well as on the coefficients Aqq, Aqq and Aqq will of order 5p. Hence, Eq. [79] implies 



A FCC/ 



that both A^^^{p,p + Sp) and A^'^"^{p,p + 5p) are of order {Spy. This is the rationale 
behind the simplifying assumption made in Sect. 15.41 

At any rate, A {p,p + 5p) and A {p,p -\- 5p) can be numerically computed 
from Eq. [75] For all values of A*' simulated in Ref. "SS], their combined effect on the 
determination of the coexistence pressure turns out to be smaller than 1% of the 
statistical error bars 1411. 



5.6 Algorithmic performance 



The simulation of the weight in Eq. H71|) requires two types of moves: single particle 
displacements, as well as changes in the volume of the simulation box. We shall use 
the short hand Elementary Monte Carlo Step (EMCS) to the combination of N con- 
secutive single-particle displacements attempts] | followed by a change attempt in the 

simulation box volume. 

The standard tools to asses the performance of a Monte Carlo algorithm [45] are 
briefly recalled in the Appendix [XI 



^* We pick at random a particle-index, say i, and try r^ — > r^ -|-(5 with 6 chosen with uniform 
probability within the sphere of radius A. We tune A to keep the acceptance above 30%. 
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Fig. 17 (Color online) Integrated autocorrelation times, defined in Eg. 1881 for Qa (top) and v 
(bottom), as a function of S (tlie Unear coordinate that labels the integration path in Fig. 1151 
where 5 = stands for the fluid minimum and 5 = 1 represents the homogeneous FCC phase) . 
Time is measured in units of EMCS. Results for all values of N equilibrated in Ref. I38| . 



Recall that we will be discussing independent simulations along the path in Fig. 1151 
We thus will be labeling them by means of a coordinate S, such that S = corresponds 
to the fluid pure phase, while 5=1 refers to the FCC minimum. 

One should like to consider the time autocorrelation functions for the components 
of the gradient field, bq^ and b(j. Yet, Eq. H73|l tells us that these correlation functions 
are identical to those of Qq{R) and C{R). Eq. (|77p suggests as well that the time 
autocorrelation function for the specific volume v is of interest. An example of these 
autocorrelation functions is shown in Fig. 1161 for the S — point. We note that v plays 
the role of the algorithmic slow mode, with a strong A'' dependence. On the other hand, 
the autocorrelation function for Qg decreases very fast, and it is barely A-dependent. 
The autocorrelation function for C is qualitatively identical to that of Qg, and will 
thus be skipped. 

The analysis is made quantitative by considering the integrated autocorrelation 
times, see Fig. 1171 We notice that the dynamics of v is considerable slower than that of 
Qq, and featureless as a function of S. Data for the specific volume scales as Tv ~ A ' 
(quite worse than standard critical slowing down in three dimensions, t ^ N ' , yet 
much better than exponential dynamic slowing-down). There is a clear anomaly in the 
behavior of r for a single simulation point in A = 2916. We briefiy comment on this 
below. 

Using these tools, it was ensured in Ref. |38] that all simulations were, at least, 
50r long. Besides, all simulations were performed twice, with different starting config- 
urations (either an ideal FCC crystal, or an ideal gas). Compatibility between the two 
sets of investigations was systematically checked [41] . 
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Fig. 18 (Color online) For systems of A^ hard spheres at their phase-coexistence pressure, we 
show, as a function of the line parameter 5, different linear combinations of the particle-density 
fluctuations, Eq. II85I I. computed for the minimal wavevectors allowed by periodic boundary 
conditions and ordered in such a way that J-"i > J-2 > J-3. For the phase-separated states all 
three J-" are of order 1 (order 1/N for homogeneous systems). The slab phase is the only one 
with J-"i — J-2 of order one. The cylinder phase is identified by J-2 — J-3 of order one. 



The anomaly at S = 0.4 for A*' — 2916 is due to the emergence of metastability. At 
this value of S we expected to find a spatially segregated state (a slab of FCC crystal 
in a liquid matrix). This state appeared indeed, but the simulation tunnels back and 
forth from it to an helicoidal crystal (in close analogy with the Monte Carlo history 
shown in the bottom-left panel of Fig. llOp . We actually performed extra simulations 
in this point, in order to determine the relative weight of each of the two metastable 
states (their eflect was carefully considered in final estimates [41j). 

These helicoidal crystals appear much more often for A'^ — 4000 and intermediate S. 
Nevertheless, selecting carefully the starting particle configuration for the simulation 
at each S, one may obtain a gradient field with a smooth 5-dependency. However, it 
is clear that these A'^ = 4000 results, although plausible, cannot be regarded as well 
equilibrated [38] . 



5.7 Interfacial free-energy 



The interfacial free energy is the free-energy cost per unit area of a liquid-to-crystal 
interface. Its computation has been rather difficult for hard spheres, different authors 
finding mutually incompatible results |90ll91ll92| . 

The tethered formalism allows as well to compute the interfacial free energy [38] . 
One should identify a point along the straight path in Fig [T5] where a slab of fluid 
surrounded by an FCC crystal is formed. Such configurations have two phase-separating 
surfaces of linear dimensions comparable with those of the simulation box. The effective 
potential difference between such a heterogeneous state and either of the two pure 
phases provides an estimate of the interfacial free energy. It turns out that a well 
formed surface appears only for N > 2048 [38]. Let us see how this can be ascertained. 
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The physical situation is as follows. When we go from the liquid to the solid, 
Fig. 1151 the homogeneous fluid becomes unstable at a value of the linear coordinate 
S oc N~ '^ ~^ ' , which means that a macroscopic droplet of crystal forms. This has 
been established for all types of first-order phase transitions '24 ,25,26,28 , and explic- 
itly verified for crystallization [38]. As S grows the mass of the crystal droplet increases, 
which costs surface energy. At a certain point, the periodic boundary conditions allow 
reducing the surface energy by turning the crystal droplet onto a crystal cylinder. At 
still larger S, the cylinder becomes a slab. Of course another three analogous geomet- 
rical transitions arise when S keeps increasing as we approach the FCC minimum. All 
six geometric transitions appeared in large enough hard-spheres systems [38J. We are 
interested in identifying systems large enough to form a slab of crystal surrounded by 
fluid. 

To follow these geometric transitions we consider the particle-density fluctuations 
quantified through 

N 2 



^('^) = ^E^'"" 



(85) 



As we are interested in the largest wavelength, we consider the smallest q allowed by 
periodic boundary conditions, ||q|| — 2tt/L, where L is the linear size of the simula- 
tion box. There are three such minimal wavevectors in a cubic box. Given a particle 
configuration, J-"i is the maximum over the three directions, J-3 is the minimum, and 
J^2 is the intermediate one. As the droplet, cylinder and slab geometries have different 
symmetries the natural order parameters are 

— Whenever the system is phase separated, (J-'i + J-2 + J-3)/3 is of order 1, (order 
1/A'' otherwise). 

— For a cylinder, two of the J-"'s are of order 1, while the J-" along the cylinder axis is 
small. Hence, J-2 — J-3 is of order 1 in the cylinder phase, but it vanishes (for large 
A'^) both in the droplet and the slab phase. 

— For a slab the only T of order 1 is that transversal to it. Hence J-"i — J-2 is of order 
1 for a slab, but not for the cylinder nor the droplet. 

All these behaviors are identified in Fig. 1181 We thus conclude that A^ > 2048 is 
sufficient to attempt a computation of the interfacial free-energy. 

To improve the numerical stability, we consider the potential difference between 
two points, S* ~ 0.5 and S — 1. For both S values, the gradient of J7jv is normal to 
the straight in Fig [T3] (actually S = 1 is a local minimum of J7jV) while S* is a local 
maximum if we restrict ourselves to the straight). Hence, the interfacial free-energy for 
the (100) crystal direction, 7100 > follows from [38] 

m^jvi/3 (^FCC-^5») ^ (86) 

A peculiarity of the tethered approach is that one may control the dependence 
of the estimate of 7100 on the actual estimate used for the coexistence pressure. One 
simply computes 7100 as a function of pressure, using Eg. 1861 as it is shown in Fig. 1191 
It turns out that the slope of the curve is of order 0.4, hence an error of order e in the 
determination of p^ results in an error of order ~ 0.4e in 7100- To our knowledge, such 
effects have not been taken into account in previous computations [90II91II92] . 
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Fig. 19 (Color online) Intcrfacial free-energy for the (100) lattice crystalline direction 7100 
as a function of pressure, for a system of N hard-spheres, for N = 2048, 2916 N = 4000. We 
estimated 7100 from Eq. 1861 



6 Conclusions 



In this work, we have explained Tethered Monte Carlo 36] , a systematic and very gen- 
eral methodology to overcome free-energy barriers in Monte Carlo simulations without 
a random-walk in the reaction- coordinate space. An admittedly weak point is that 
a deep physical insight is needed in order to identify a suitable reaction coordinate. 
Nevertheless, it is worth mentioning that this problem may be bypassed using real 
replicas [42| . 

The power of the method is illustrated in two famous problems of modern theo- 
retical physics, the diluted antiferromagnet in an external field (DAFF), the physical 
realization of the random field Ising model, and the crystallization of hard spheres. In 
both cases, the tethered approach allows to go beyond the state of the art. 

Let us stress that Tethered Monte Carlo is not limited to reproducing canonical 
mean values, but it is also suitable to address different physical questions that one 
would never ask in a canonical setting due to overwhelming technical difficulties. Our 
study of the pseudo phase coexistence in the DAFF, see Sect. 14.71 is one example. 
Furthermore, for disordered systems, tethered Monte Carlo provides a redefinition of 
the quenched disorder, less vulnerable to (trivial) violations of self- averaging [691139) . 

We remark as well that, from the algorithmic point of view, we have been rather 
unsophisticated. We used local Metropolis moves that (in the case of the DAFF) were 
complemented with parallel tempering. Clearly enough there is much room for im- 
provements on this side (see e.g. [37] for an implementation of a cluster algorithm). 

We also note that, in our implementation, both for the DAFF and for hard spheres 
there is a maximum system size beyond which one cannot equilibrate within reason- 
able computer time. Clearly enough new free-energy barriers appear, which are not 
adequately described by the chosen reaction coordinates. These could be treated by 
extending the number (or the type) of tethered quantities. 
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A Assessing thermalization in Monte Carlo simulations 

This appendix is intended as a brief and handy reference on thermalization in Monte Carlo 
simulations, providing some definitions that are used throughout the paper (see |45l for a de- 
tailed discussion). Most of the definitions are standard, but subsection lA.ll describes a powerful 
thermalization criterion that may be new even for the Monte Carlo specialist. 

As a general rule, the thermalization of a Monte Carlo simulation should be assessed 
through the temporal autocorrelation functions 1451 . If 0{t) is the measured value of the 
observable O at time t during the simulation, then 

Co(i,m) = ([O(0)-(O>^,][O(i)-(O>^])^,, poit,m) = ^^l^'"!] . (87) 

™ Co(0, m) 

From this function we obtain the integrated and exponential autocorrelation times 

rint,o{rh) = - + J2po{t,rh), (88) 

^ t=i 

rcxp,o(m) = lim sup — , (89) 

t^oc -\og\po{t,m)\ 

Toxp(m) = supTcxp,o("i)- (90) 

o 

In general, the autocorrelation function can be expressed as a sum of exponential^^ 

po(i,A) = ^A,e-*/-». (91) 

i 

The exponential time, then, is the largest of the t^ and characterizes the time needed to equili- 
brate a certain observable (or the whole system). The integrated time indicates the minimum 
time difference so that two measurements of some observable O can be considered indepen- 
dent (i.e., uncorrelated) . Notice that if the decomposition II91II contains a single exponential, 
f"exp,0 = f"int,Oi This is a useful observation because the exponential time is harder to measure 
accurately than the integrated one. 



^^ Eq. II91I I holds only for an aperiodic dynamics that fuUfils detailed balance. Violation of 
either of these two condicions will result in the appearance of complex eigentimesTi . If the 
dynamics is aperiodic, the corresponding eigenvalue of the generator of the Markov process, 
e~^''^', will be smaller than one in absolute value. Hence for these modes, the correlations 
could take the form of a damped oscillation. A well known example is the autocorrelation 
function for the magnetization in the Wolff single-cluster dynamics for the Ising model (see, 
e.g.. Fig. 3-5 in |44l). Indeed, cluster algorithms verify balance, but not detailed balance, 
because the corresponding generator of the Markov chain is formed by alternating steps of two 
detailed-balance fuUfiling dynamics, namely bond-updating and spin-updating. In the case of 
a detailed-balance dynamics, the eigenvalues e"^'"^* are real but they could still be negative. 
Having said so, in most applications (and certainly in the ones considered here), oscilating 
behavior is not observed within statistical errors. 



42 

A.l Thermalization in parallel tempering simulations 

The computation of the tq requires a much longer simulation time than what is needed to 
thermalize the system. This is not a problem for ordered systems, since extending the simula- 
tion reduces the statistical errors. In disordered systems, however, the errors are dominated by 
sample-to-sample fluctuations, so an efficient CPU time allocation should run each sample for 
a time long enough to thermalize it, but no longer. Therefore, the computation of autocorrela- 
tion times in these cases has traditionally been abandoned as impractical and thermalization 
is typically assessed by less rigorous methods such as the time evolution of disorder-averaged 
observables. 

The use of parallel tempering, however, provides a way to compute the autocorrelation 
times in relatively short simulations. The key is analyzing quantitatively the temperature ran- 
dom walk of each participating configuration §3, . Reference 2Q includes a detailed discussion 
of this method, here we shall only give a brief sketch of it. 

Let us consider a parallel tempering simulation with Nt participating replicas of the 
system. We define fi{t) as the function indicating the temperature index of replica i at time 
t. By 'temperature index' we mean the temperature mapped onto {1,2, ..., Nt}, where 1 
corresponds to the lowest temperature and Nt to the highest. We can now compute the 
autocorrelation function II87I I of each fi (notice that {fi) = {Nt + l)/2). The Nt correlation 
functions thus obtained can be averaged, so they collectively have the precision of a much 
longer simulation. 
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